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The previous investigation on Rayleigh-Benard convection of a dilute classical gas [T. Kita: J. 
Phys. Soc. Jpn. 75 (2006) 124005] is extended to calculate entropy change of the convective tran- 
sition with the rigid boundaries. We obtain results qualitatively similar to those of the stress-free 
boundaries. Above the critical Rayleigh number, the roll convection is realized among possible 
steady states with periodic structures, carrying the highest entropy as a function of macroscopic 
mechanical variables. 
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I. INTRODUCTION 

In a preceding paper,— we performed a statistical me- 
chanical investigation on Rayleigh-Benard convection of 
a dilute classical gas based on the Boltzmann equation. 
We specifically calculated entropy change through the 
convective transition for the case of the stress- free bound- 
aries as a function of macroscopic mechanical variables. 
We thereby tested the validity of the principle of maxi- 
mum entropy proposed for nonequilibrium steady states? 2 * 
The present paper extends the consideration to a more 
realistic case of the rigid boundaries^ 



II. FORMULATION AND NUMERICAL 
PROCEDURES 



Our starting point is Eq. (46) of ref. [l|: 
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where V • j'^ 1 - 5 ) =0 and the units are described in §3.1. 
We adopt the condition of the rigid boundaries along z, 
i.e., j (1 - 5) =0 at z = ±l/2. Combined with V-j( 1 ' 5 )=0, 
it yields the boundary conditions: 
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the system, i.e., the total particle number, energy, and 
energy flux along z. These conditions lead to eq. (42) of 
ref.Q], i.e., 
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where ATh c = 2jq ^ /innP^ denotes the temperature dif- 
ference between z — ±1/2 that would be realized in the 
heat-conducting state. 

We solve the above equations with the method devel- 
oped by Peschi^ First, the boundary conditions of eq. ([2]) 
is treated with the Galerkin method, 5 i.e., by expanding 
every z dependence in terms of some basis functions sat- 
isfying the boundary conditions. Specifically, the basis 
functions for eq. (|2a|) are obtained from the second-order 
differential equation S" = -AS with 5(±l/2)=0 as 

S n {z) = \/2smmr(z + l/2) (n = 1, 2, • • • ). (4a) 
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They satisfy {S n \S n ') = J_ 1/2 S n (z)S n >(z)dz = 5 nn >. On 

the other hand, those for eq. (|2b|) are constructed from 
the fourth-order differential equation — k 4 C with 
C(±l/2) = C"(±l/2) = 0as 
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where fc„(>0) is determined by 
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with jf' 5 " 1 denoting the xy components. As for the hor- 
izontal directions, we consider the region — L/2 <x,y < 
L/2 with L^> 1 and impose the periodic boundary con- 
ditions. We also fix macroscopic mechanical variables of 
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Note k n « (n+l/2)7r for n^$> 1. The quantity A n is the 
normalization constant: 
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so that (C„|C„') =<5 rm '. The functions {C„(z)} may be 
called Chandrasekhar functions^ 

Now that appropriate basis functions are obtained, we 
expand and j' 15 ^ in eq. U) as 
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Here fc^ = component is excluded in the expan- 
sion of j'^ 1 - 5 ) to seek only periodic current distribu- 
tions in the xy plane. We have also incorporated into 
eq. (|5a[) the fact that the temperature is uniform at 
z = ±1/2. Now, V • ji^ 1 - 5 ) = is transformed into 
ik± ■ j±(k±,n)+J2 nf {Sn\C! n ,)~j z (k±,n') = 0. It hence 
follows that j±(k±, n) can be expressed generally as 
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On the other hand, eq. ([3]) is transformed into 
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Let us substitute eq. with eq. © into eq. (fT]) and per- 
form space integrations using the orthonormality of the 
basis functions. We thereby obtain algebraic equations 
for the expansion coefficients as 
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with (f\g)^j%dxj^dyj^ /2 dzr(r)g(r). 

Finally, entropy characteristic of convection is obtained 
by substituting eq. (|5a|) into eq. (50) of ref. [l| as 
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with Ti and AT given by eq. ©. 

It follows from the stability analysis for the heat- 
conducting stated that the critical Rayleigh number R c 
is determined from eq. {SJ by setting the nonlinear terms 
and time derivatives equal to zero. The relevant insta- 
bility originates from the linear coupled equations for 
T(k,n) and j z (k,n). Eliminating T(k,n) in favor of 
j z {k, n), we obtain the equation for R c as 



detA = 0, 
where matrix A is defined by 
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The critical Rayleigh number R c corresponds to the min- 
imum value of in eq. (|10[) as a function of k±. 
Equation (|10|) is solved by approximating A by a finite 
dimension of n c xn c , and the convergence is checked by in- 
creasing n c . Choosing n c — 4 already yields an excellent 
result of R c = 1.708 with kj_ = 3.116 = fc c £ 

The nonlinear terms become relevant in eq. ([5]) for 
R<--V > # c . They are evaluated for given expansion 
coefficients T(k±,n) and j(k±,n) as follows. We first 
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construct TW(r) and j( 15 )(r) by eqs. The fast 

Fourier transform (FFT)& is used in this procedure to ob- 
tain the xy dependence. We then perform the space dif- 
ferentiations numerically in the xy plane and analytically 
along the z direction. We specifically use the following 
formulas of 0(h 6 ) in the xy plane: 

d x f(x,y)n ]T ^[f(x+3ah,y)-9f(x+2ah,y) 

CT=± 

+45f(x+ah,y)}, (13a) 

d 2 J(x, y) « { E [2/(^ + 3^, y)-27f(x+2ah, y) 

+270/(x+afc, y)} - 490/(s, y) J. (13b) 

The quantity d xy f{x, y) is obtained with eq. (|13a|) by av- 
eraging the derivatives performed in different order. The 
nonlinear overlap integrals are evaluated finally, where we 
again use the FFT in the xy plane. On the other hand, all 
the calculations along z are performed by preparing the 
relevant overlap integrals in advance, e.g., (S n \S n >S n ») 
and {S n \C n i S' n ,,) for eq. (|8a[) . and performing the sum- 
mations over n! and n". 

Time evolutions of eq. (J5J) are calculated as follows. 
We first multiply eqs. ((HEJ) and (JHcJ) by (OT 1 ) n "n and 
[fci + (n^) 2 ]- 1 with {0) nn > = k 2 ± 6 nn ,-{C n \C^), respec- 
tively, and perform summation over n for eq. (|8bp ■ Time 
integrations are then carried out numerically by treating 
df/dt = g as /(t+At)« f(t)+g(t)At. A disadvantage 
of this simple method is that we have to make At small 
enough to avoid an explosion in the numerical time in- 
tegration. One may alternatively use the split-step inte- 
gration scheme developed by Pesch which approximates 
df/dt = Lf+g as / (t±At)«e^ At /(t) + e^ 2 [3g(t) - 
g(t — At)] At/2, thereby treating the linear part Lf ex- 
actly. This latter scheme removes the explosion at the 
expense of larger numerical errors to make a rapid time 
integration possible; the extra computational time for di- 
agonalizing L in the calculation of e— At is negligible in 
the whole numerical procedures. 

We here focus on periodic solutions of eq. ([1]) in the 
xy plane and express rj_ = siai + S202, where a\ = 
{ai x , ai y , 0) and a 2 = (0,a2,0) denote the basic vec- 
tors. Accordingly, we adopt the periodic boundary con- 
dition for the region spanned by A/jOi and ^20,2 with Afj 
(j = 1,2) a large integer. The above theoretical frame- 
work can also be used in this case with a minor modifi- 
cation. Indeed, we only have to perform the change of 
variables (x,y) — > (sx,S2) in the xy integrations of the 
nonlinear terms. Those integrations have to be carried 
out now only over the unit cell of < s\ , S2 < 1 ■ The cor- 
responding wave vector k± is given by k± = t\b\ +^2^2, 
where b\ = 2ix{a,2 X e z )j\(a\ x 02) • e z ], 62 = 2ir(e z x 
ai)/[(cii x 02) -e z ], and £j denotes an integer. The linear 
stability analysis for the heat-conducting state suggests 
that the stable solution satisfies ~ ]&2| ~A; C = 3.116. 



The parameters in eq. ([8]) are chosen the same as those 
used for the free boundaries, i.e., eq. (37a) of ref.[0, which 
correspond to Ar at 273K under atmospheric pressure. 
We also fix the heat-flux density jg at z = — 1/2 so that 
the temperature difference AT\ IC = IK is realized between 
z = ±1/2 in the heat-conducting state. The Rayleigh 
number is controlled by changing the thickness d. 

Practical calculations of eq. ([8]) are performed as fol- 
lows: We first multiply eqs. (gll)-® by 10 s , 10 10 
and 10 10 , respectively, and rewrite them in terms of 
f'{kj_,n) = 10 3 f(fc ± ,n) and j'(k±,n) = 10 5 j(kj_,n) to 
obtain equations of 0(1). The summations over n are 
truncated at a finite value n c , whereas Appx discrete 
points are used to perform FFT for each direction in 
the xy plane. As for periodic structures, we investigate 
the three candidates: the roll, the square lattice and the 
hexagonal lattice with |6i| = | £>2 1 ~ k c . We then trace 
time evolutions of the expansion coefficients until they 
all acquire constant values. Choosing At < 0.005, n c >4 
and iVpFT ^ 2 4 yields excellent convergence for the cal- 
culations presented below even with the simplest time- 
integration scheme. The initial state is chosen as the 
conducting state with small fluctuations T'(k±, 1) ~ 10~ 2 
for the basic harmonics k±. The constants AT and Ti 
are updated at each time step by using eq. ([7]) . Also eval- 
uated at each time step is entropy measured with respect 
to the heat-conducting state: 

AS = - S£ ] , (14) 

where is given by eq. © and S$ = -5(AT hc ) 2 /48. 
We thereby trace time evolution of AS simultaneously. 
The above procedure is carried out for each fixed periodic 
structure. 

One of the advantages of the present approach is that 
we only have to change the basis functions to study 
other boundary conditions. For example, the case of 
the stress-free boundaries can also be treated within the 
present framework by simply changing C n (z) — > S n (z) 
and S n (z)— s-%/2cosn7r(z + l/2) in the expansions of eqs. 
([5b|) and (|5c)) . respectively. This replacement has been 
checked to reproduce the results obtained in ref. [I] appro- 
priately. 

III. RESULTS 

We now present numerical results on the rigid bound- 
aries, which turn out to be qualitatively the same as those 
of the stress-free boundaries^ 

Figure [Tj shows time evolution of A5* for the Rayleigh 
number = 1.2i? c . The letters r, s and h de- 

note (r) roll, (s) square and (h) hexagonal, respectively, 
distinguishing initial conditions; they are exactly the 
same as those for the stress-free boundaries.— Writing 
T'(k±,n) = T'[lx,l2, n] and introducing the angle 9 by 
9 = cos~ 1 (bi • b2), those initial conditions are given 
explicitly as follows: (r) f '[±1,0,1] = 1.00 x 10~ 2 /%/8 
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FIG. 1: Time evolution of entropy measured with respect 
to the heat-conducting state for = 1.2i? c . The letters 

r, s and h denote roll, square and hexagonal, respectively, 
distinguishing initial fluctuations around the heat-conducting 
solution; see text for details. The final state of t > 50 is 
the roll convection, whereas the intermediate plateaus of s 
and h correspond to the square and hexagonal convections, 
respectively. 
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FIG. 3: Profile of the average temperature variation T(z) in 
the roll convection normalized by the temperature difference 
ATkc in the heat-conducting state. The Rayleigh numbers 
are R { ~ 1) =R C , 1.2R C , 2.0R C , 5.0R C and 10.0i? c from top to 
bottom on the left part. 
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FIG. 4: The length |bi|/fc c of the stable roll convection as a 
function of the normalized Rayleigh number R^ 1 ' /R c . 



FIG. 2: Time evolution of entropy AS. The four curves cor- 
respond to the different Rayleigh numbers: i?' -1 - 1 = 1.2i? c , 
2.0.Rc, 5.0i? c and 10.0J? C . The initial state is the heat- 
conducting state with the fluctuation T'[±1,0, f] = 1.0 x 
10~ 2 /\/8 and |6i| = fe c , whereas all the final states are the 
roll convection. The broken line near the top indicates the 
upper bound of AS. 



with |6i| = k c ; (s) f'[±l,0, 1] = 1.01 x 1(T 2 /V8 an d 
f'[0,±l,l] = 0.99 x W- 2 /V8 with |6i| = |6 2 | = k c and 
9 = tt/2; (h) f'[±l, 0, 1] = T"[0, ±1, 1] = 1.00 x lQ~ 2 /\/8 
and f'[l,l,l] = f '[-1,-1,1] = 1.01 x 10- 2 /V8 with 
\bi | = 1 62 1 = k c and 9 = 2ir/3. We observe clearly 
that entropy increases monotonically in all the three 



cases to reach a common final value of the roll convec- 
tion. The intermediate plateaus seen in s and h corre- 
spond to the metastable square and hexagonal lattices, 
respectively, with (s) T' [±1,0,1] - T"[0,±1,1] and (h) 
f'[±l, 0, l]~f"[0, ±1, l]~f"[±l, ±1, 1]. 

Figure [2] displays time evolution of AS for four dif- 
ferent Rayleigh numbers, all developing from the initial 
fluctuation f'[±l, 0, 1] = 1.00 x 10" 2 /v% with |&i|=Jfe c . 
Each final state is the roll convection. Compared with 
the case of the stress-free boundaries^ we observe an en- 
hanced oscillatory behavior for R = 5.0i? c and 10.0i? c 
after the first rapid increase of AS*. 

Figure [3] shows profile of the average temperature vari- 
ation T(z) along z in the roll convection for five different 
Rayleigh numbers. Figure 3] plots |i>i|/fc c corresponding 
to the maximum of AS as a function of R^/ R c . A gam 
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the basic features are qualitatively the same as those of 
the stress-free boundaries^ 

Thus, we have seen that the principle of maximum en- 
tropy proposed in ref. is satisfied through the Rayleigh- 
Benard convective transition of a dilute classical gas even 
in the realistic case of the rigid boundaries. 
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